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ABSTRACT 

This thesis considers two models for the computation of 
position finding confidence, one of which utilizes a bivariate 
normal distribution of region, and the other a chi-square 
distribution. The two models are based on different 
assumptions, these are explained and explored. A computer 
simulation model is presented which utilizes both position 
finding models under varying conditions. 
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I. 



INTRODUCTION 



Bearings are taken on an emitter (target) of unknown 
geographic location by n stations of known location. The 
best point estimate (BPE) of the target is then calculated, 
along with a confidence, or search, region within which the 
target may be said to be located with a given probability. 

This constitutes the classic radio direction finding position 
location problem. 

This classic problem is part of a general class of 
position estimation problems whose basic approach is identical, 
which is the calculation of a point location from a number of 
lines of position. 

Variant problems in this class occur when lines of position 
are established from ranges, or relative times of arrival of 
an emission, or when the object to be located takes the 
measurements on the reference stations. 

While this paper examines the classic problem, much of the 
methodology applies to this more general class of problems. 

The BPE is the primary result of a position location 
solution. The confidence region amplifies this information 
in two ways. The size of the confidence region serves as a 
measure of the reliability of the BPE, a larger region 
implying a more tenuous location. The region also describes 
an area in which the target will be located with specified 
confidence. The geometric details of this region provide the 
more likely locations for search. 
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II. COMPUTATION OF CONFIDENCE REGIONS 



A typical [6] sequence of operations for a position 
estimation by computer is as follows: A first point esti- 

mate is made and a plane surface tangent to the earth is 
established at that point. Remaining computations are made 
in the plane rather than the surface of the earth. Unusable, 
or wild, bearings are rejected and a BPE made utilizing the 
remaining bearings. A confidence region is computed using 
the same bearings. 

More than one mathematical method is available for most 
of the steps [3, 4, 8, 9, 11] . The two methods used herein 
for confidence region computation are designated the 
bivariate normal (BVN) and the chi-square (x 2 ) • 

The Bivariate normal region [3, 4] is the most commonly 
used. This region is the family of concentric ellipses 
ax 2 - 2bxy + cy 2 = -2 log (1 - P) 
where P = probability that the target is found within the 
ellipse 



x, y = coordinates in local (to the target) reference 
system 

£ cos 2 6 i 
a = E J_ 

j=l o? 

n sin 0 • cos 0 . 

b = Z 2 1 

j=l a\ 



n sin 2 



c = l 



j=l ot 
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0j = bearing observed by station j 

a ? = variance of displacement of line of bearing 

n = number of bearings 

The bearings from the stations are assumed to be normally 
distributed with a mean of zero, and are displaced parallel 
to the true bearing at the site of the BPE . The confidence 
region lies in a plane tangent to the surface of the earth 
at the BPE as calculated by the least squares method [3, 4]. 

The BVN confidence region requires several assumptions: 

(1) The position lines are great circles on the 
surface of the earth. 

(2) Measurement error is Gaussian with mean zero. 

(3) Displacement of position lines from the true 
position line is parallel to the true position 
line . 

(4) The earth is flat in the vicinity of the 
confidence region. 

These last two assumptions represent approximations that 
are made to simplify the calculation of the confidence region. 
Thus, the calculated region is an approximation. The parallel 
bearing displacement is an approximation of the actual dis- 
placement of a point on the bearing by a bearing error. From 
this, the probability distribution of the points is developed 
[4]. Placing the resultant curves in a plane adds an 
additional distortion, which grows with the size of the 
confidence region [7]. 
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The chi-square confidence region [3, 10] is the family 



of concentric curves 



X 



2 _ 



n 

z _L 

1=1 S2 



< 0 . - 






v;here 

p 

Sj = bearing variance of station j 
0_. = bearing observed by station j 
3 = true bearing from station j to a point 

j 

(x,y) in the curve 
n = number of bearings 

p 

X has a chi square distribution with n degrees of 
freedom. Of the above listed assumptions required by the 

p 

BVN method, only (1) and (2) are required for the x region. 

The region lies on the surface of the earth centered at the 
2 

minimum X point, rather than in a plane. 

Both methods require the estimated variance of the 
bearings to be used in the calculation. The effect of an 
incorrect estimation could be much larger than the errors 
of approximation. Estimates in practice are a function of a 
stations past performance .plus the direction finder operator's 
evaluation of the reliability of each bearing he observes and 
reports . 

The assumption of Gaussian bearing distribution eliminates 
the case of wild, or outlying bearings which occur in practice, 
largely through operator error. As a Gaussian distribution 
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is basic to both the BVN and x methods, elimination of wild 
bearings from the computation is necessary, though the theory 
of each method would be affected by their assumed existence. 

Wild bearing rejection methods are available [6, 7]. The use 
of such methods is usually limited to four or more bearing 
situations. With two bearings, no basis for detection of 
wild bearings exists. For most three bearing situations, when 
the existence of at least one wild bearing is indicated, it 
is impossible to identify it, or them. With four bearings, 
one wild bearing can be seen to miss the BPE area formed by 
the other three . 

In practice [6, 7] the Outlier detection and rejection 
method is used, but a Gaussian distribution is still assumed 
for those remaining. 

Daniels and Vajda [3] object to the use of x regions on 
the grounds that it may occur that no points can be found to 
satisfy the resultant equations. 

Read [10] examines the simplifications required in 
calculating the BVN regions, and compares the two methods. 

The x^ regions are expected to be too small, and the BVN too 
large, though the full effect of the simplifications by Kukes 
and Starik [4] is not known. 

Comparison of the two methods by computer plotting 

techniques shows considerable difference in size. Additionally, 

the effect of the approximations used in the BVN method is 

apparent. While BVN region is elliptical in all cases, the 
2 

X region may assume an egg shape. 
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A barrier to the use of x regions is the lack of a way 
of reporting them. A BVN region can be described by its 
center, angular orientation, and length of semi axes. No 
such simplifying parameters exist for x regions. One 
possible approach would be to attempt a conversion into a 
partially equivalent rectangle, square, or circle as is often 
done for BVN regions [1, 6, 7], Since these regions are 
always larger than the computed region, contain areas of 
lesser probability, and omit areas of greater probability, 
any advantage of using one method over the other would 
probably be lost thereby. 

A thorough investigation of the two methods is required, 
in order to determine the actual probabilities represented 
by each type of confidence region. The differences indicate 
that at least one of the methods has error. A computer 
simulation was written to compare them. 
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III. COMPUTER SIMULATION MODEL 



The appended computer program provides a means for 
simulation of the two confidence region methods for identical 
input data. 

Inputs to the program are station latitude, longitude, 
and standard deviation of bearing distribution; target 
latitude 'and longitude. 

The true bearings and distance for each station and 
target are computed and a bearing error generated and added. 
This portion contains parts of the Bergh program in Pope [8] 
and is a program used by LT Paul Winkler , modified here to 
include more general positions. 

A first-point estimate of the position is made by Pope's 
[9] vector method, and the final fix by the least squares 
method of Kukes and Starik [4] in a plane tangent to the 
earth at the first-point estimate. This portion was used by 
Lunde [5] but contained errors since corrected. 

A bivariate normal confidence region is computed by the 
method of Kukes and Starik, and orientation angle and 
length of major and minor simi axes generated. Additionally, 
output by graphic plot is available. 

A confidence region is then computed and a graphic 

2 

plot generated. No parameteric output for the x region is 
available, as no method yet exists for describing the shape 
of the region. This portion was originally written by LT 
Winkler, and since modified to include more general positions. 
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The subroutines GAMA, for the normalized incomplete gamma 
function, and MTRMAP , the graphic plot, are from the library 
of the W. R. Church Computer Center at the Naval Postgraduate 
School . 

The graphic plots are on a grid measured in degrees of a 
great circle, resembling a Mercator projection at the equator. 
Hence, visual distortion of the regions occurs as their 
location approaches the poles, but the geographical points 
of the x regions are correct. The geographical points of 
the BVN regions are correct, except for the distortion caused 
by the flat earth assumption. 
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IV. SUGGESTED SIMULATION PROCEDURE 



An undetected error in previous work invalidated all 
data generated. Time constraints prevented regeneration of 
the data after correction of the error. 

The behavior of the regions generated by the two methods 
differs with the number of stations in the net, configuration 
of the net, standard deviation of the error used, and relative 
location of target with respect to the net. Thus, a general 
case may not exist. A valid comparison would include 
sufficient combinations to detect all possible patterns of 
behavior . 

A useful measure of effectiveness for the regions would 
be the difference in the probability level in each of the two 
regions occupied by the true position of the target. 

The standard deviation used in generating bearing error 
is also used in computing BPE and confidence regions. A 
minor modification to the program would permit the examination 
of the consequences of incorrectly estimating the actual 
bearing distribution of a station, by computing these with 
a different standard deviation that that used to generate 
the error. 

The errors generated are from a normal distribution, and 
include no wild bearing. A means of detection for wild 
bearings is not included in the program. Such a routine 
may be added to the program, along with a generation routine 
for wild bearings. 
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The modifications above would permit examination of more 
realistic situations through the use of operational data, 
along with a model of the actual net which supplied the data. 
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V. AN ALTERNATE APPROACH 



An underlying assumption in the methods discussed is that 
there is a uniform prior probability for the target location. 
That is, the mathematical solution of the BPE and confidence 
regions are inputs to the solution of the operational 
question of target location. Prior knowledge of probable 
target location does not affect the solution. This approach 
is intuitively correct generally. 

Operationally, there are two uses of position finding 
information : 

(1) To determine the location of a target of known 
or unknown identity; 

(2) To identify a target from a set of known targets 
by its location. 

This second use introduces prior information into the 
problem, and occurs usually when attempting to identify a 
fixed land-based station. 

Butterly [2] describes a method for incorporation of 
information of this type in producing a Bayesian position 
estimate, though he implies that its use could be more general. 

A further modification of the program would permit a test 
of the efficacy of this method by comparing the number of 
correct identifications when using the Bayesian method with 
the number when selecting the target in the lowest level of 
the confidence region. Both the BVN and x methods could be 
used in each. 



14 



non 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 



THIS FORTRAN PROGRAM CALCULATES RADIO DIRECTION FINDING POSITION 
ESTIMATES AMD CONFIDENCE REGIONS. SEVERAL OPTIONS IN OUTPUT ARE 
AVAILABLE. AND ARE DOCUMENTED BY COMMENT CARDS IN THE PROGRAM. 



INPUTS TO THE ^RAM q ARE : STATIQNS (0N g CARD, FORMAT 12) 

LONGITUDE, AND STANDARD DEVIATION Or EACH STATION 
EACH, FORMAT 3F10.4) 

TARGETS (ONE CARD, FORMAT 12) 

AND LONGITUDE OF EACH TARGET (ONE CARD EACH, 



NUMBER OF 
LATITUDE, 
(ONE CARD 
NUMBER OF 
LATITUDE 



FORMAT 2F10 .4- ) 

5 - SIZE OF GRAPHICAL PLOTS IN 



DEGREES (ONE CARD, FORMAT F10.4) 



LIMITS - STATIONS 8 
TARGETS 8 



SIGN CONVENTION: 



LATITUDES 
LONGITUDES 
BEARI NGS 



NORTH 0 TO +9 0.0 

SOUTH 0 TO -90.0 

EAST 0 TO -180.0 

WEST 0 TO +180.0 

EAST OF NORTH 0 TO +180.0 

WEST OF NORTH 0 TO -180.0 



PART ONE - INITIALIZATION AND INPUT 



DIMENSION BFTA( 10, 10) 

DIMENSION STEEST (8,8) 

DIMENSION STALAT(IO) 

DIMENSION TARLON ( 20) 

WORM ( 13, 13) 

PJ(1 1) , Z E T A ( 1 0 ) 

B 1 ( 8 ) , B 2 ( 8 ) , v 

At ( 8 ) , BE ( 8 ) , CE ( 8 ) , DE ( 8 ) 
Y ( 13 ,13) 

SDRADt 8) 

C (13 ,13) 



DIMENSION 
DIMENSION 
DIMENSION 
DIME NS ION 
DIMENSION 
DIMENSION 
DIMENSION 



GAMMA (8,8) , „ 

STALON(IO), STDEV(IO), T ARLAT ( 20 ) 
THE TA ( 10,20), DIST( 10,20) 



EE( 8) , F E ( 8 ) 



REAL*4 T ( 24 ) 

RADDEG=. 017453293 
I X=4672 1 
PI=3. 141592 
RADI US = 10800. O/PI 
EX = 0.0 



READ STATION POINTS 
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non oooooooooooo noon noon 



READ ( 5 i 100) NS T A 
100 FORMAT { 12 ) 

2 READ ( 5 * lol^STALATlI ) * STALCNll), STDEV(I) 
102 FORMAT (3F10. 4) 

READ TARGET POINTS 
READ ( 5 » 100 ) NT AR 

DO 1 I = 1 *NT AR , , 

1 READ ( 5 » 102) TARLAT(I)* TARLON(I) 

READ GRAPH SIZE 

READ ( 5 * 102 ) CGRAPH 



PART TWO - GENERATION OF BEARINGS AND BEARING VARIANCES 

THIS SECTION OF PART TWO COMPUTES BEARING ANGLES AND DISTANCES 
USING NAPIER'S ANALOGIES 

DO 10 1 = 1* NSTA 
DO 1C J = 1 * NT AR 
AA=AB S { STALCN! I ) — T ARLON { J ) ) 

IF( AA .GT .130 .0 ) AA = 360.0-AA 
GAMMA! I * J) =AA 

CHECK FOR SAME MERIDIAN CASE 

IF(AA.GT.O.Ol) GO TO 105 
BB=0 . 0 

I F (STALAT ( I ) .GT .TARLAT { J ) ) BB= 180.0 

d'iST^ I ] J )=ABS( STALAT! I ) -TARLAT! J ) ) *RAD I U S*RADDEG 
GO TO 23 

105 CONTINUE 

CC=ABS! AA-180. 0) 



16 



ooooooo noon on o 



I F ( CC . GT. 0. 001 ) GC TO 106 
DD=0.0 

SUMLAT=STALAT( I ) +TARLAT ( J) 

I F ( SUMLAT. LT.O. 0) DD=180.0 
THET A( I , J)=DD 
EED=180 . 0-ABS { SUMLAT) 
DIST(I»J)=EED*RADI US*RADDEG 
GO TO 23 
106 CONTINUE 



CONTINUE SOLUTION OF TRIANGLE 



W VALUE = GAMMA ( I » J )*. 5*RADDEG 
COT 1 = COS (WVALUE ) / S IN ( WV ALU E ) 
BRAVO = 90.0 - TARLAT(J) 

ABLE = 90.0 - STALAT(I) 



COS1 

C0S2 

SIN1 

SIN2 



COS( (ABLE 
COS( (ABLE 
S IN ( (ABLE 
SI N( (ABLE 



- BRAVO)*. 5*RADDEG) 
+ BRAVO )*. 5*RADDEG) 

- BRAVO)*. 5*RADDEG) 
+ BRAVO)*. 5*RADDEG) 



Y A=COT l*COS 1 
YB=COT 1* S I N 1 



ATA=ATAN2(YA,CGS2) 
ATB=ATAN2( YB »S I N2) 
B3B=I ATA-ATBJ/RADDEG 



S8B=ANGLE » NOW COMPUTE PROPER SIGN 



SL=STAL0N( I ) 

TL=TARLON(J ) 

IF(SL.LT.O.O) SL=360.0+SL 
I F ( TL. LT.O. 0) TL=369.0+TL 
TL =TL-SL 

IF(TL.LT.O.O) TL=T L+360 . 0 
IFtTL.LT. 180.0) BBB=-BBB 
THET A { I . J ) =BBB 

THE TA=TRUE BEARING FROM STATION TO TARGET 



DISTANCE FROM STATION TO TARGET USING LAW OF COSINES 
GAMMA ( I , J) = GAMMA ( I » J ) * RADDEG 
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TARLAK J) =TARLAT { J ) *RADDEG 

TEMP A i l SI N (TARLAT ( J T ST AL AT ( I ) )+COS(TARLAT(J) J^COSISTALATI I ) ) 

l*COS(GAMMA( I » J ) ) 

DFDIST = ARCOS(TEMP) 

DIST(I.J) = DFDISKRADIUS 

DIST = DIST ANCE IN NAUTICAL MILES 

TARLAK J)=TARL AT ( J) /RADDEG 
STALAT ( I )=STAL AT ( I ) /RADDEG 

23 CONTINUE 



Tut r c p r t t nf ' 1 (IF PART TWO ADJUSTS THE STANDARD DEVIATION TO BE USt-D 
1W RFARING e’p^OR COMPUTATION FOR THE DISTANCE BETWEEN STATION AND 
TARGET. THE SOURCE OF THIS COMPUTATION IS NOT KNOWN. 

SD = STDEV ( I ) 

R = DIST (I , J )/ 100.0 
IF ( R— 54. 0) 20,20,22 t 
22 SD = SD* (R/1108.0 - R)) 

GO TO 32 

2 J IF (R-10.0) 24,26,26 

24 IF (R-4.C) 28, 30,30 
28 SD = P 1/2 .0 

30 SD =°SD*( .0204*R*R - . 402*R + 3.0) 

26 SD IsDK ,O0071*R*R - .3213*R + 1.1598) 

32 STEEST { I » J ) = SD*SD*RADDEG*RADDEG 
10 CONTINUE 

253 FORMAT 6 !/// 136X.45H BEARING ANGLE FROM ITH STATION TO JTH TARGET,/ 

*WRIT E( 6 , 256 ) ( ( I,J,THETA(I,J),J=1,NTAR) , I =1 , NSTA) 

256 FORMAT ( 39X , 13HF ROM STATION ,12, 

111H TO TARGET ,I2,4H — ,F10.4) 

255 FORMAT 6 (//fi OX, 4 OH DISTANCES FROM ITH STATION TO JTH TARGET,//) 

WRITE(6,256)((I,J,DIST(I,J),J=1,NTAR),I-1,NSTA) 
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C THIS SECTION OF PART TWO COMPUTES BEARING ERROR FROM A SIMULATED 
C NORMAL DISTRIBUTION USING A RANDOM NUMBER GENERATOR , AND ADDS THE 
C ERROR TO THE TRUE BEARING. 

C 

WRITEI6, 260) 

260 FORMAT I // ,50X , 1 5HN0RMAL VARIATES,//) 

C 

DO 40 1=1, NSTA 
DO 40 J = 1 , NT AR 

STEEST ( I , J ) = STEESTd , J) /RADDEG 
CALL GASS( IX, STEESTd ,J) , EX, X) 

BETA(I , J ) =THET A ( I , J ) +X 
C 

WRITEI6, 256)1, J, X 
C 

40 CONTINUE 



OUTPUT ANGLES WITH ASSOCIATED ERROR 
WRITE (6 ,200 ) 

200 FORMAT ( // 6X , 6H STA LAT , 8X, 6HSTAL0N, 8X , 6 HTARLAT , 8X.6HTARLON, 
18X, lOHTHETAt I, J ) , 4X, 5HSTDEV) 

WRITE ( 6,259) ( ( STALAT (I ) , S T A L ON ( I ) , T ARLAT ( J) , T ARLON ( J ) , 

1 BET A ( I , J ) , STE E ST( I , J ) , J = 1,NTAR), 1=1, NSTA) 

259 FORMAT (/,6 (F12.5.2X ) ) 



PART THREE - COMPUTATION OF POSITION ESTIMATE, OR FIX 



DO 3321 J1=1,NTAR 

WR I T E ( 6 , 220 ) T ARLAT ( J 1 ) , T ARLON ( J 1 ) 

220 FORMAT ( 71H 1 THE FOLLOWING BEARINGS WERE TAKEN ON A TARGET LOCA 

1 T ED AT LATITUDE, F12.5,11H LONG I TUDE , F 12 . 5 , // ) 

WRITE (6,300) 

300 FORMAT ( 10X, 8HLATITUDE , 12X , 9H LONG I TUDE , 1 4X , 5HTHET A , 1 1 X , 1 2 HT RU E BEAR 
1ING, 10X, 7HSTD DEV,/) 

WR I T E I 6 , 22 1 ) ( STALATd ) ,STALON( I ) , BET A { I , J 1 ) , T HET A ( I , J 1 ) , 

1 ST DEV ( 1 ) , 1 = 1, NSTA ) 

221 F0RMAT(5(10X,F10.5) ) 
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THIS SECTION OF PART THREE CALCULATES A FIRST POINT ESTIMATE OF 
THE FIX USING POPE'S VECTOR METHOD (1971) 



COMPUTE BEARING PLANE INTERSECTIONS 
WRITE (6, 701) 

701 FORMAT ( / / ,44X , 26HBEARI NG LINE INTERSECTIONS. 
1/.40X, 8HSTAT IONS. 4X.8HL ATI TUDE ,4X, 
29HL0NGITUDE.// ) 

DO 88 I =1. NSTA 
C 

B1 ( I ) = STALAT (I ) * RADDEG 
B2 ( I ) = STALON ( I ) * RADDEG 
C 

SLAT = S I N ( B 1 ( I ) ) 

CLAT = C OS ( B1 ( I) ) 

SLON =- S I N ( B2( I ) ) 

CLON = COS (B2( I ) ) 

SBNG = SIN(BETA(I,J1)*RADDEG) 

CBNG = COS(BETA( I ,J1)*RADDEG) 

C 

XE = CLAT * CLON 
YE = CLAT * SLON 
ZE = SLAT 
C 

AE ( I ) = CBNG- SLON - SBNG* SLAT*C LON 
BE f I ) = -CBNG*CLON - SBNG* SLAT -SLON 
CE ( I ) = SBNG*C LA T 
C 

D E ( I ) = BE ( I )*ZE - C E ( I ) *YE 
E E ( I ) = CE { I ) *XE - A E ( I ) * ZE 
F E ( I ) = AE(I)*YE - BE ( I ) *XE 
C 

88 CONTINUE 
C 

USUM =0.0 
VSUM =0.0 
WSUM = 0.0 

c 

DO 89 I =2. NSTA 
ILIMIT = 1-1 
DO 90 J=1 , ILIMIT 
C 

UE = BE(I)*CE(J) - CEII)*BE(J) 

VE = C E ( I ) *AE( J ) - AE ( I ) *C E ( J ) 

WE = A E ( I ) * B E ( J ) - BE ( I ) * AE ( J ) 
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c 



T I 
TJ 



UE*DE( I ) 
U E* DE ( J ) 



+ VE*EE( I) + WE-FE(I) 
+ VE*EE ( J ) + WE*FE( J ) 



C 



IF !TI*TJ .LT. 0.0) GO TO 90 
IF til .GT. 0.0) GO TO 91 



UE = -UE 
VE = -VE 
WE = -WE 



COMPUTE CENTROID OF INTERSECTIONS 

91 USUM = USUM + UE 
VSUM = VSUM + VE 
WSUM = W SUM + WE 

UMAG=SQRT ( U E*U E+VE*VE +WE*W E ) 

ULAT=ARSIN( WE/UMAG) /R ADD EG 
ULON=-AT AN2 ( VE * U E ) /RADDEG 

WRITE! 6,702) I , J,ULAT,U LON 
702 FORMAT (40X, 12, 5H AND , I 2 , 3X, F 8. 3 , 4X , F8. 3 ) 

ULAT, ULON = BEARING PLANE INTERSECTIONS 

90 CONTINUE 
89 CONTINUE 



CONVERT CENTROID TO LATITUDE AND LONGITUDE OF FIRST POINT 

XI LON =-ATAN2( VSUM, USUM) /RADDEG , wrs 

X1LAT= A RSIN ( W SUM/SQRT { USUM* USUM+VSU M* VSUM+W SUM* W SUM ) )/R 



ESTIMATE 

ADDEG 



222 

223 



WRITE!6,222 ) 

WRITE! 6,223) X1LAT, 
FORMAT!//, 1GX , 24HT H E 
FORMAT ( 1 OX , 8 HX 1 L AT = 



XI LON 

FIRST POINT ESTIMATE,//) 

, F10.5 . 5X» 8HX1L0N = ,F10.5,//) 



THIS SECTION OF PART THREE TRANSFERS THE PROBLEM TO A PLANE TAN^cNT 
TO THE EARTH AT THE FIRST POINT ESTIMATE, THEN COMPUTES A POSITION 
ESTIMATE USING THE LEAST SQUARES METHOD OF KUKES AND STARIK. 



ASUM=0 . 0 
BSUM=0 . 0 
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CSUM=0 . 0 
DSDM=0. 0 
ESUM=0 .0 



DO 8999 L = 1 . NS TA 
SDRADl L)=STDEV ( L )*RADDEG 
8999 CONTINUE 

DO 9002 I = 1 » NS TA 
XSI=6ETA(I«J1)*RADDEG 
XS I =AB S I XS I ) 

XS I =STAT ION BEARING IN RADIANS 



COMPUTE BEARING IN LOCAL REFERENCE PLANE USING 
AND L AW OF SINES 

RHO= ( STALON ( I ) -X 1L0N }*RADDEG 

RH0=A3S I RHO ) # niJO 

IF(RHO.GT.PI) RHO=( 2. O^PI ) -RHO 

CHECK FOR DUE NORTH OR SOUTH BEARING 

IFtRHO.GT. 0.001) GO TO 345 

IF( STALATC I ) .LT.X1 LAT) ZETA{ I ) =0.0 
GO TO 347 

345 CONTINUE 
RR=ABS( RHO— P I ) 

IF (RR.GT .0 .001 ) GO TO 346 
ZE T A ( I ) =P I 

POLE = ST A LAT ( I ) +X1LAT 
I F ( POL E . LT .0 .0 ) ZET A ( I ) = 0 • 0 
GO TO 347 

346 CONTINUE 

CONTINUE SOLUTION OF TRIANGLE 

C E E= ( 90 .-STALAT ( I ) )*RADDEG 
DENOM=. 5*( XSI *-RHG) 

RUM11= .5*(XSI-RH0) 

RUM12=. 5*CEE 

ANUM=S IN { RUM 11 )«TAN(RUM12) 

DENA=SI N(DENOM) 

ARKTN1=ATAN2( ANUM.DEMA) 

RT A=ABS ( DENOM-P I / 2 ) 



NAPIER' S ANALOGI ES 
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RTB=ABS(RUMll-PI/2) _ 

IF(RTA.GT.C.OOl) GO TO 355 
BNUM=CO S ( RUM11 ) ’•‘TAN! RUM 12) 

DENB=0 .001 

GO TO 357 , „„ 

355 IFIRTB.GT. 0.001) GO TO 356 
BNUM=0 .001 

DENB=COS ( DENOM ) 

GO TO 357 

356 CONTINUE , 

BNUM=C0S(RUM11 ) *TAN( RUM12) 

DENB=C0S ( DENCM ) 

357 CONTINUE _ , 

ARKTN2=ATAN2(BNUM,DtNB) 
AAEE=ARKTN1+ARKTN2 
IFIXSI.LT. 0.001) XSI=0.001 

CDD=SIN(CEE)*SIN( XSI ) „ , 

I FtAAEE.LT .0 .001 ) AAEE=0 .00 1 
EEE=SIN(AAEE) 

FFF=DDO/EEE ^ 

IF( FFF.GE.l .0) FFF = 0. 999999 
CCEE=AR SI N ( FFF ) 

ZET A( I ) =P I-CCE E , 

IFIAAEE.LT. CEE) ZETA(1)=CCEE 7CTA , t > 
I F ( B ET A ( I » J 1 ) .LT.0.0) ZETA( I ) =-ZETA( I ) 



347 CONTINUE 
ZET A= LOCAL I ZED BEARING 



TAKE NEGATIVE OF ZETA TO CONVERT TO NORMAL TRIG 
ZE TA ( I ) =— Z E T A ( I ) 



CALCULATE PERPENDICULAR DISTANCE FROM LOCALIZED 
CENTER OF COORDINATE SYSTEM 



YPRIM=(90.— XILAT )*RADDtG 

AAEE=ABS( AAEE) 

Y D I ST= (YPRIM-AAcE) 

PJ( I )=YDI ST*SI N( ZETAI I ) ) 



SUM TERMS FOR FOLLOWING CALCULATIONS 



SIGN CONVENTION 
BEARING TO 
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ASUM=ASUM+ (COS ( ZETA(I))*COS(ZETA(I)))/ ( S DRAD ( I )*S DRAD( I ) ) 
BSUM=BSUM+(S1N (ZETA( I ) )^'C05(ZETA( I ) ) ) / ( SDP.AD ( I ) * SDR AD ( I ) ) 
CSUM=CSUM+ ( S I N ( ZET A ( I ) ) 7; S IN ( ZETA( I) ) )/ (SDRAD( 1 }*SDRAD( I) ) 



C 

9002 CONTINUE 
C 

DSUM2DSUMiPJ(Kli((BSUM*SIN(ZETA(K) i-CSUM*COS(ZETA(K)))/(SDRAD(K)* 

1 ESUM=E SUM+P J ( K ) *( ( A SUM* SIN(ZETA(K} ) -BSUM*C0S ( ZET A ( K) ) )/ (SDRAD(K)* 
1SDRAD(K) ) ) 

C 

9006 CONTINUE 



CALCULATE LEAST SQUARES ESTIMATES AND CONVERT TO DEGREES 

C0EF=1 .0/ ( (ASUM*CSUM)-(BSUM*BSUM) ) 

XO=COEF*DSUM 
YO=COEF*ESUM 
S0K0=-2.*AL0G( . 1 ) 

XFST=X1LQN+X0/RADDEG 
YEST=X1L AT+YO/ RADDEG 

XEST, YEST = LEAST SQUARES ESTIMATE OF POSITION 



WRITE ( 6 ,9007) YEST, XEST 
5007 FORMAT (• ' , ‘THE LEAST SQUARES 

1* DEGREES LATITUDE' //43X , F7 , 



ESTIMATE OF 
2,' DEGREES 



POSITI ON IS ' ,F6. 2 , 
LONGITUDE'// / ) 



PART FOUR - COMPUTATION OF CONFIDENCE REGIONS 



THIS SECTION OF PART FOUR CALCULATES A BIVARIATE NORMAL CONFIDENCE 
REGION (DANIELS, KUKES AND STAR1K) 

CALCULATE ANGLE BETWEEN MAJOR AXIS OF ELLIPSE AND MERIDIAN 

AA1=2.*BSUM 
AA2=CSUM— ASUM 
TW0ANG=ATAN2(AA1,AA2) 
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BANG=TW0ANG/2. 0 

IF(BANG)651,651,652 

651 BANG=BANG+l PI/2. 0) 
GO TO 653 

652 BANG=BANC-( PI/2.0) 

653 CONT INUF 

ANG= BANG/RAD DEG 



CALCULATE SEMIAXES FOR ELLIPSE OF .9 PROBABILITY 

t 

IFIBANG.GT. 1.56) BANG=1.56 
IFIBANG.LT. -1.56) BANG=-1 . 56 

SDENOM = CSUM + BSUM*T AN ( BANG) 

FALSE=ABS(SDENOM) 

SAMA JO= SORT! SO KO /FALSE ) 

SDEN0= ASUM-BSUM^TAN ! BANG ) 

F ATS =ABS ( S DENO ) 

SAMINO=SQRT (SOKO/FATS ) 

SAMINO = S/,MINO/ RA CD EG 
SAMAJO=SAMA JO/RADDEG 



WRITE! 6, 9081) ANG , SAMAJO , SAMI NO 
9081 FORMAT! • 10X ANGLE BETWEEN MAJOR AXIS AND MERIDIAN 

134X, ’SEMI-MAJOR AXIS = * , F6 .2 , / / , 34X * 'SEMI-MI NOR AXIS 



F 6 . 2 / / 

F 6 . 2 / / / ) 



THIS SECTION CF PART FOUR PRINTS A GRAPHIC PLOT OF A BIVARIATE 
NORMAL CONFIDENCE REGION. THE SIZE IS DETERMINED BY THE INPUT 
•CGRAPH'. USE THE SMALLEST VALUE CF CGRAPH THAT WILL DISPLAY 
THE REGION, TO MINIMIZE DISTORTION. 

THIS SECTION MAY BE SKIPPED IF THIS OUTPUT IS NOT DESIRED. JUMP 
FROM THIS POINT TO THE BEGINNING OF THE CHI-SQUARE SECTION, BELOW . 

LOAD PROBABILITY ARRAY FOR BIVARIATE NORMAL 

XOR=XE ST 
YOR = YE S T 
DO 9011 J=l, 13 
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EXVAL={ C GRAPH/2. 0) - ( { ( J-l) / 1 2. 0 ) *CGRA PH ) 

EXVAL=EXVAL*RADDEG 
DO 9012 1=1,13 

YVAL=(CGRAPH/2.0)-( ( { I- 1 ) / 12 . 0 ) *CGRAPH ) 

YVAL=YVAL*RADDEG 

H00=BSUM*E XVAL^YVAL— ( A SUM/ 2. ) * EXVAL* EX VAL- ( CSUM/2 .0 J *YV AL*YV AL 
PR0B= 1 . 0— EXP ( H00 ) 

WORM { I , J) =PROB 
9012 CONTINUE 
9011 CONTINUE 
C 

N= 13 
M=13 
AZ=0 . 0 
BZ=0 . 0 
BND=0 . 0 
AMI N = 0. 0 
I JT= 0. 0 
I CON= 1 
I GR = 1 

DO 9080 LL = 1 , 2 4 
9080 T ( L L ) = 0 . 0 



PRINT THE CONFIDENCE REGION PLOT 

CALL MTRMAP ( WORM , N, M, T , BND , AZ , BZ , AM I N , I J T , ICON, IGR } 
WRITE(6,591) 

591 FORMAT {/// ,45X, ' BIVARIATE NORMAL CONFIDENCE REGION') 

W R I T E ( 6 , 509) CGRAPH ,CGR APH 

509 FORMAT t///,20X,37H THE SIZE OF THE CONFIDENCE PLOT IS ,F10.4, 
114H DEG. LAT. BY ,F10.4, 12H DEG. LONG. ,/,40X, 

241 H THE FIX POINT IS LOCATED A1 (007, 007). ) 



THIS SECTION OF PART FOUR CALCULATES A CHI-SQUARE CONFIDENCE 
REGION, AND PRINTS A GRAPHICAL PLOT. NO OTHER OUTPUT IS AVAILABLE. 
ADDITIONALLY, LEVELS OF CONSTANT VALUE OF Q MAY BE PLOTTED (DANIELS) 
THIS SECTION MAY BE SKIPPEO BY JUMPING FROM THIS POINT TO THE END OF 

EITHER THE Q OR CHI-SQUARE PLOTS MAY BE SKIPPED 
JUMPING AROUND THE APPROPRIATE 'CALL MTRMAP...' CARD 
THIS SECTION. 



THE MAIN PROGRAM. 
SEPARATELY BY 
AT THE END OF 



X1LAT=YEST 
XI LON=XE ST 
X1LAT = X1LAT 
X1L0N = X 1 LON 



+ CGRAPH/2.0 
+ CGRAPH/2.0 
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DO 501 K = 1,13 
DO 502 L = 1,13 
V ( K , U ) = 0.0 
C 

DO 571 1=1, NST A 
G=ABS ( S TALON ( I ) -XI LON) 

I Ft G.GT .180 .0 ) G=360.0— G 
IF (G.GT. 0.001) GO TO 610 
ANGLE= 0.0 

IF(STALAT( I) .GT.X1LAT) ANGLE= 180.0 
GO TO 620 
610 CONTINUE 

ABLE=9 0.0— STALAT (I ) 

BRA V0=9 0. 0-X1 L AT 

WVALUE=G*. 5*RADDEG 

COT1 = COS ( WVALUE J/SIN (WVALUE) 

COS1 = COSt ( ASLE-BRAVO)*.5*RADDEG) 

COS 2 = COSt (ABLE+BRAVC)*.5*RADDEG} 

SIN1 = SIN( (ABLE— BRAVO)*. 5-RADDEG) 

SIN2 = SINt (ABLE+BRAVO)*.5*RADDEG) 

RR = COT 1*S IN 1 
SS=AT AN2 { RR , S I N2 ) 

TT=COT 1*C0S1 

UU = AT AN2 (TT ,C0S2 ) 

ANGLE=( UU-SSJ/RADDEG 
620 CONTINUE 

SL=STA LON ( I ) 

T L = X 1 L C N 

IFtSL.LT .0.0) SL=360.0+SL 
IF(TL.LT.O.J) TL=360 .0+TL 
TL=TL-SL 

IFtTL.LT .0.0 ) TL=T L+360 . 0 
IFtTL.LT.lt, 0.0 ) ANGLE=-ANGLE 
CHIDIF=ANGLE-BETA U ,J1 ) 

Y ( K , L } = Y { K , L ) + CHIDI F*CHIDI F/STDEVt I )/STDEV( I ) 
P AR 1=N S TA/2 . 

PAR2=Y( K,L )/2. 

CALL GAMA ( FAR1 , P AR2 , G AM , B , 0 .0 ) 

C ( K, L ) = 1 .—GAM 
571 CONTINUE 
C 

502 X1LCN s X1LCN - CGRAPH/12.0 

X1LCN = XI LON + CGRAPH- 13. 0/12.0 
501 X1LAT = X1LAT - CGRAPH/12.0 



SET UP THE CONTOUR SUBROUTINE 
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N = 13 
M = 13 
AZ = 0.0 
BZ = 0.0 
BND = 0.0 
AMIN = 0.0 
IJT = 0 
ICON = 1 
I GR = 1 

DO 504 II = 1,24 
504 T ( I I ) = 0.0 

PLOT THE 0 LEVELS 

CALL MTRMAPtY, N,M,T, BND, AZ,BZ, AMIN, IJT, ICON ,IGR) 
WRITE (6, 592) 

592 FORMAT (///,45X, 'LEVELS OF CONSTANT VALUE OF Q* ) 
WRITE! 6,509) CGRA PH ,CG RAPH 

PLOT THE CHI-SQUARE REGION 

CALL MT RMAP ( C , N , M , T , BND , AZ , 3Z , A MI N , I JT , I C ON , I GR) 
WRITE(6,593) 

593 FORMAT! /// ,45X CHI -SQUARE CONFIDENCE REGION') 

K R I T E ( 6, 5 09) CGRAPH ,CGRAPH 

C 

3321 CONTINUE 

C 

STOP 

END 



SUBROUTINE GAS S ! I X , S , AM , V ) 

A = 0.0 

DO 50 1=1,12 

CALL WRANDU ( IX , IY, Y ) 

IX = IY 
50 A = A+Y 

V = ( A-6 .0 ) -S+AM 

RETURN 

END 
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SUBROUTINE WRANDUt I X , I Y »Y) 

IY = IX*65539 

IF { IY) 5,6.6 _ , 

5 IY = IY + 2147483647+1 

6 Y = IY 

Y = Y* . 4656613 E— 9 

RETURN 

END 
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